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1.  Introduction  and  Background 


The  Soldier  Protection  Sciences  Branch  of  the  US  Army  Research  Laboratory 
(ARL)  relies  heavily  on  numerical  models  of  the  brain  and  surrounding  tissues  to 
understand  load  distributions,  which  will  allow  the  design  of  efficient  and  effective 
protective  equipment  for  the  Army.  A  recent  study  of  interest  to  the  ARL 
community  used  both  2-  and  3-dimensional  (2-D)  and  (3-D)  models  to  simulate  the 
pressures  and  velocities  within  the  brain  when  it  was  cyclically  rotated  about  the 
axis  of  the  neck  at  a  targeted  rate.'  The  models,  meshed  using  Eulerian  (fixed-grid) 
particle-based  methods  in  the  code  Uintah,  were  validated  via  an  experiment  in 
which  human  subjects  laid  down,  placed  their  heads  in  a  device,  and  rotated  their 
heads  back  and  forth;  tagged  MRI  was  used  to  collect  relevant  data.  Although  the 
magnitudes  and  propagation  of  velocities  between  the  empirical  data  and  simulated 
results  in  the  study  were  similar,  the  predicted  pressure  distribution  in  the  brain  was 
shifted  in  the  results  of  the  human  experiments.  Noticing  that  the  complicated 
interface  between  the  brain,  cerebrospinal  fluid  (CSF),  and  skull  was  not  well 
resolved  in  the  study  (properties  were  averaged  among  the  3  materials  along  the 
boundary),  we  hypothesized  that  the  discretization  scheme  may  skew  such  results. 
Our  objective  in  this  study  is  to  quantify  the  effects  that  discretization  methodology 
has  on  such  model  results  by  creating  and  analyzing  a  2-D  brain  model  with 
simplified  geometry  under  3  different  discretization  schemes  and  thus  under  3 
different  methods  of  resolving  the  brain/CSF/skull  interface.  We  conducted  our 
study  using  both  rotational  boundary  conditions  and  a  planar  blast-loading 
boundary  condition  for  each  discretization  case. 

2.  Methods 


To  create  models  of  the  skull,  cerebrospinal  fluid,  and  brain,  we  relied  on  the 
software  ALE3D  from  Lawrence  Livermore  National  Laboratories.^  ALE3D 
possesses  an  internal  mesh  generator  that  allows  finite  element  meshes  to  be  created 
and  simulations  run  in  Lagrangian,  Eulerian,  Arbitrary  Lagrangian/Eulerian  (ALE), 
and  embedded  Lagrangian/Eulerian  schemes.  In  a  Lagrangian  simulation,  the  mesh 
is  tied  to  the  material  as  it  moves  or  deforms;  Eulerian  simulations  keep  the  mesh 
fixed  while  allowing  material  to  advect  through  it;  ALE  methods  move  the  mesh 
independently  from  the  material  motion,  often  depending  on  how  distorted  the 
mesh  becomes  over  time;  and  the  newly  developed  embedded  method 
(implemented  via  the  ALE3D  FEusion  package)  allows  coupling  of  fully 
Lagrangian  and  fiilly  Eulerian  objects  within  the  same  simulation. 
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We  simplified  the  geometry  of  our  brain  model  by  reducing  its  shape  to  a  2-D 
ellipse  with  a  semi-major  axis  of  approximately  103  mm  and  a  semi-minor  axis  of 
83  mm,  surrounded  by  a  7-mm-thick  skull.  For  the  Lagrangian  model,  the  CSF  was 
accounted  for  using  a  frictionless  laminar  slide  surface  between  the  brain  and  skull. 
The  Eulerian  model  contained  a  1-mm-thick  CSF  layer  that  was  incorporated  via 
mixed  material  zones  in  the  mesh;  in  mixed  material  zones,  material  properties  are 
averaged  by  volume  in  each  element  that  contains  multiple  materials,  and  the 
boundaries  between  materials  are  reconstructed  after  each  time  step  based  on 
volume  fractions.  The  embedded  model  comprises  a  Lagrangian  skull  and  a 
Eulerian  brain,  and  the  CSF  was  represented  via  a  frictionless  slide  surface  between 
the  two.  In  the  pure  Lagrangian  and  pure  Eulerian  rotational  models,  a  void 
background  was  used  due  to  the  added  simplicity  it  gave  to  the  problem;  in  all  other 
cases,  a  Eulerian  mesh  of  air  was  used  as  the  background.  The  Lagrangian  blast 
loading  problem  required  that  embedded  methodology  be  used  to  resolve  the 
boundary  between  the  air  and  the  head.  Figure  1  gives  examples  of  each 
discretization  type. 
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Fig.  1  Example  meshes  for  each  discretization  case.  The  meshes  shown  use  Lagrangian, 
Eulerian,  and  embedded  methods  from  top  to  bottom. 
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The  material  properties  for  the  air,  CSF,  brain,  and  skull  were  obtained  from  Dr 
Yolin  Huang.  ^  The  air  was  modeled  as  a  gamma  law  gas  with  coefficient 
(;k  -  1)  =  0.4.  The  brain  was  modeled  as  an  incompressible  Mooney-Rivlin  material 
with  the  follow  equation  of  state  for  Cauchy  stress: 

a  =  f  [o)  (o)(B-i)’]  =f[co  B'-(l-  o))((B2)'-  tr  (B)  B)  ] .  (1) 

The  skull  model  used  was  a  linear  elastic  model  with  a  shear  modulus  of 2,660  MPa 
and  a  bulk  modulus  of  3,891  MPa.  The  units  used  throughout  our  simulations  for 
material  properties,  coordinates,  and  results  were  based  on  grams  (g),  millimeters 
(mm),  milliseconds  (ms),  and  megapascals  (MPa).  These  are  different  than  the  units 
traditionally  used  with  ALE3D  for  blast  simulations. 


The  head  meshes  in  the  rotational  models  were  rotated  using  an  angular  velocity 
load  curve  applied  to  a  node-set  on  the  skull.  In  the  Lagrangian  model,  the  node¬ 
set  chosen  was  the  very  outside  surface  of  the  skull.  In  the  embedded  simulation, 
all  sets  of  nodes  inside  the  skull  except  for  the  inner  and  outer  surfaces  of  the  skull 
were  used.  The  Eulerian  model  used  a  filled  circular  mesh  object  that  surrounded 
the  skull  and  was  made  of  the  same  material  as  the  skull  to  rotate  the  head  by 
applying  the  angular  velocity  load  curve  to  the  set  of  all  nodes  within  a  circular 
radius  outside  of  the  skull  but  within  the  circular  object  (Fig.  2).  This  method  of 
rotation  was  used  because  the  Eulerian  method  prevented  rotation  of  a  noncircular 
node-set. 


DB:  shaire  1  _eHipseEuler_  1 2_CI64.00000 
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Fig.  2  Loading  node-set  for  Eulerian  rotational  problem.  The  dark  shaded  area  around  the 
skull  is  the  area  to  which  the  angular  velocity  loading  was  applied. 
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To  evaluate  our  rotational  models,  we  ran  each  simulation  for  520.0  ms  with  a 
loading  curve  represented  by  Fig.  3.  In  this  time,  the  head  rotates  60°  to  the  left, 
120°  to  the  right,  and  60°  left  back  to  the  center  position.  The  rate  was  selected 
based  upon  android  research  performed  by  Robert  Fitzpatrick,"^  which  suggested 
that  the  maximum  rate  of  safe  rotation  of  a  human  head  is  467°/s.  For  the  embedded 
grid  simulations,  we  were  only  able  to  run  the  model  for  40.0  ms  as  opposed  to 
520.0  ms  because  of  the  long  run  time  of  the  simulations;  so  the  loading  curve, 
although  the  same  shape  as  the  one  for  the  other  2  models,  produced  much  higher 
velocities.  Because  of  the  time  constraints  of  the  project,  no  attempt  was  made  to 
optimize  the  run  time  of  the  simulations.  The  3-dimensionality  of  the  embedded 
simulations  and  the  Lagrangian  blast  simulation  did  not  appear  to  significantly 
affect  run  time. 
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Fig.  3  Load  curve  for  rotational  models 

The  force  applied  in  the  blast  loading  models  was  simulated  by  applying  velocity, 
density,  and  energy  inflow  conditions  to  the  air  node-set  on  the  extreme  left  edge 
of  the  model,  340  mm  to  the  left  of  the  outer  left-most  edge  of  the  skull.  The  inflow 
conditions  were  obtained  by  collecting  time  history  values  from  a  modified  version 
of  Dr  Joel  Stewart’s  ShockTube  model,  which  simulated  a  TNT  explosion  initiated 
at  the  apex  of  a  rigid  right-circular  cone  that  was  14  m  long.^  The  time  histories 
were  collected  2  m  radially  inward  from  the  widest  part  of  the  cone  and 
approximately  halfway  between  the  top  edge  of  the  cone  and  the  midpoint  axis  of 
the  cone  in  the  vertical  direction.  The  ShockTube  simulation  ran  for  5.0  ms,  and 
time  histories  were  collected  throughout  the  simulation.  The  time  history  data  used 
were  then  restricted  to  the  time  in  which  a  blast  was  actively  moving  past  our 
selected  tracer  particle  location. 

Our  blast  loading  simulations  were  thus  run  for  2.5  ms.  Figure  4  depicts  the  inflow 
conditions  used. 
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Velocity  Load  Curve,  Blast  Simulations 


Energy  Load  Curve,  Blast  Simulations 


Time,  ms 


Density  Load  Curve,  Blast  Simulations 
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Fig.  4  Velocity,  energy,  and  density  inflow  conditions  for  planar  blast  loading  problem 

Because  of  the  planar  symmetry  of  the  blast  loading  problems,  the  problem  was  cut 
in  half  along  the  head’s  major  axis,  and  only  the  upper  half  of  the  problem  was  run 
(Fig.  5). 
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Fig.  5  Example  of  planar  blast  loading  problem.  The  dimensions  of  the  air  in  the  simnlation 
are  900  mm  in  the  x-direction  by  450  mm  in  the  y-direction.  The  example  shown  is  pnre 
Enlerian. 

Data  collection  was  performed  using  6  tracer  particles  located  at  different  locations 
inside  the  brain  in  each  mesh;  the  locations  of  the  tracer  particles  are  shown  in  Fig. 
6.  For  the  purposes  of  our  simulations,  the  right  side  of  the  brain  is  the  side  of  the 
brain  closest  to  the  right  side  of  the  screen  when  viewed  in  a  postprocessor  during 
the  first  time  cycle.  The  top  side  of  the  brain  is  the  side  closest  to  the  top  of  the 
screen  when  viewed  in  a  postprocessor  during  the  first  time  cycle. 

In  one  of  the  Lagrangian  rotational  simulations,  2  additional  time  history  tracers 
were  used  to  calculate  the  positions  over  time  of  a  pair  of  adjacent  master  and  slave 
surface  nodes  on  the  slide  surface  to  determine  how  much  slippage  of  the  brain 
occurs  during  rotation. 
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Fig.  6  Locations  of  tracer  particles.  The  particles  are  located  5, 20,  and  50  mm  in  from  the 
right  and  top  inner  edges  of  the  skull.  The  example  shown  is  a  Lagrangian  rotational  model. 
The  red  and  green  materials  represent  the  hrain  and  skull,  respectively. 

Each  type  of  simulation  was  ran  at  several  different  mesh  resolutions;  the  average 
element  sizes  in  the  head  in  each  case  are  given  in  the  following  Table.  The  number 
of  elements  given  for  each  simulation  is  the  total  number  of  elements,  including 
background  void  or  air  elements,  in  that  simulation.  Average  element  sizes  are 
calculated  for  the  rotation  simulations  and  the  Eulerian  and  embedded  blast 
simulations  by  finding  the  percentage  area  of  the  head  within  the  mesh  and  dividing 
the  area  of  the  head  by  the  total  number  of  mesh  elements  multiplied  by  this 
percentage.  Average  element  sizes  for  the  Lagrangian  blast  simulation  were 
determined  by  finding  the  exact  number  of  elements  inside  the  head  via  temporarily 
removing  the  background  air  and  dividing  the  area  occupied  by  the  head  by  the 
number  of  elements;  this  was  done  since  the  mesh  resolution  of  the  skull  and  brain 
in  the  Lagrangian  blast  simulation  was  much  higher  than  the  mesh  resolution  of  the 
air.  Since  the  embedded  simulations  and  the  Lagrangian  blast  simulation  had  to  be 
represented  in  quasi-3-D  with  a  0.3-mm  third-dimensional  thickness,  the  element 
sizes  given  are  the  hypothetical  sizes  of  their  2-D  counterparts. 
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Table.  Average  element  sizes  within  the  head  by  simulation 


3.  Rotational  Model  Analysis 


3.1  Results  across  Discretization  Methods 

The  distributions  of  the  pressures,  von  Mises  equivalent  stresses,  and  effective 
strains  in  the  brain  at  fixed  angles  of  rotation  were  drastically  different  among  the 
discretization  schemes,  as  shown  in  Figs.  7-9.  Again,  the  embedded  simulation  was 
run  on  a  smaller  time  scale  than  the  Lagrangian  and  Eulerian  models.  Also,  the 
embedded  model  was  rotated  through  atmospheric  pressure  to  satisfy  the 
requirement  of  the  embedded  calculations;  thus,  the  range  of  the  pressures  and  Von 
Mises  equivalent  stress  for  the  embedded  simulation  shown  in  Figs.  7-9  has  been 
adjusted  for  an  initial  internal  pressure  of  0.101325  MPa.  The  effective  strain  rate 
was  given  by  the  following  equation: 

eff  strain  =J^(£ii^+  £22^+2812^)  •  (2) 
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Fig.  7  Effect  of  discretization  method  on  pressure  distribution.  The  rotation  angles  shown 
are  10.1°  from  equilibrium,  0.9°  back  toward  equilibrium  from  the  inflection  point 
(displaced  59.1°  from  equilibrium),  and  44.3°  back  toward  equilibrium  from  the  inflection 
point  (displaced  15.7°  from  equilibrium). 
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Fig.  8  Effect  of  discretization  method  on  von  Mises  equivalent  stress  distribution.  Like  in 
Fig.  7,  the  rotation  angles  shown  are  10.1°  from  equilibrium,  0.9°  back  toward  equilibrium 
from  the  inflection  point  (displaced  59.1°  from  equilibrium),  and  44.3°  back  toward 
equilibrium  from  the  inflection  point  (displaced  15.7°  from  equilibrium). 
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Fig.  9  Effect  of  discretization  method  on  effective  strain  distribntion.  Like  in  Figs.  7  and  8, 
the  rotation  angles  shown  are  10.1°  from  eqnilibrium,  0.9°  back  toward  equilibrinm  from 
the  inflection  point  (displaced  59.1°  from  equilibrium),  and  44.3°  back  toward  equilibrium 
from  the  inflection  point  (displaced  15.7°  from  equilibrium). 

One  immediately  visible  result  is  that  the  pressure,  stress,  and  strain  waves 
propagating  through  the  brain  take  on  different  shapes  and  frequencies  in  each  case. 
In  the  Eulerian  and  embedded  grid  models,  the  pressure  waves.  Fig.  7,  cause  peak 
pressures  to  arise  more  interior  to  the  brain  than  in  the  Lagrangian  case.  The 
embedded  grid  model,  although  unfortunately  under-resolved  because  of  the  time 
constraints  of  this  study,  gives  the  closest  results  to  the  Uintah  fixed-grid 
simulation.  The  Lagrangian  model,  oddly,  appears  not  to  have  a  clearly  defined 
pressure  wave  but  has  pressures  and  stresses  that  gradually  build  up  over  time,  as 
would  be  expected  in  a  problem  involving  rigid  body  rotation.  Figure  10  plots  the 
distance  between  nodes  in  a  master/slave  node  pair  over  time;  the  distance  is 
coherent  with  what  would  be  expected  of  a  fluid  rotating  within  a  solid  shell. 
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Distance  between  Master  and  Slave  Nodes  over 
Time,  44,460  Mesh  Elements 
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Fig.  10  Distance  between  master  and  slave  nodes  over  time 

An  analysis  of  what  causes  pressures  to  build  instead  of  reflect  in  the  Lagrangian 
model  is  ongoing,  but  the  preliminary  investigation  has  shown  that  the  pressure 
increase  can  be  completely  attributed  to  volume  changes  in  the  brain.  The  change 
in  energy  at  all  time  steps  in  our  simulations  is  zero,  so  none  of  the  volume  changes 
are  due  to  thermal  expansion.  It  may  be  the  case  that  the  volume  changes  are  due 
to  roundoff  error  in  our  boundary  conditions  that  causes  the  skull  to  change  shape. 

Another  immediately  visible  result  is  that  the  pressures  in  the  brain  in  the 
Lagrangian  model  are  at  least  an  order  of  magnitude  higher  than  the  pressures  in 
the  Eulerian  model.  The  magnitudes  of  pressures  observed  in  the  embedded  grid 
model  cannot  be  subject  to  comparison  with  the  other  2  methodologies  because  of 
the  shorter  timescale  over  which  rotation  was  applied.  The  scattered  pressures 
within  the  skull  in  the  Eulerian  model  arise  from  the  exterior  skull  material  used  to 
rotate  the  head. 

One  notable  result  is  that  the  highest  von  Mises  stresses  and  effective  strains  in  all 
cases  occur  very  near  to  the  center  of  the  brain.  This  differs  from  the  experimental 
results  and  brings  up  a  number  of  questions  as  to  why  this  phenomenon  arises  in 
our  models.  The  most  likely  explanation  for  this  is  the  geometry  chosen  for  our 
simulations.  The  simplified,  elliptical  geometry  chosen  for  our  simulations  neglects 
the  folds  and  lobe  patterns  within  a  real  brain  structure.  These  folds  and  lobes 
provide  a  mechanism  similar  to  hundreds  of  slide  surfaces  via  which  shear  stresses 
arising  from  the  rotation  of  the  brain/skull  interface  are  prevented  from  propagating 
far  into  the  center  brain.  Without  those  folds  and  lobes  providing  a  built-in 
mechanism  for  slip,  the  greatest  shearing  will  occur  at  the  center  of  the  brain. 
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A  concerning  phenomenon  present  in  the  embedded  grid  model  that  is  not  well 
depicted  by  Fig.  7  is  slight  pressure  loss  and  subsequent  recovery  over  the  course 
of  the  simulation.  Visual  analysis  of  results  from  the  simulation  show  that  by  the 
time  the  brain  reaches  an  inflection  point,  the  average  pressure  in  the  brain  has 
dropped  by  at  least  0.02  MPa  (the  pressure  has  decreased  to  at  least  80%  of  initial 
pressure).  All  pressure  is  regained  by  the  end  of  the  simulation;  however,  the 
distortion  in  results  that  this  generated  requires  an  understanding  of  what  in  the 
embedded  grid  method  is  causing  this  pressure  flux.  One  hypothesis  is  that  the 
reconstruction  of  the  brain/skull  boundary  in  the  embedded  calculations  is  subject 
to  roundoff  errors  when  rebuilding  element  volumes.  Small  losses  in  skull  volume 
or  stretching  of  the  skull  in  the  radial  direction  (and  thus  a  miniscule  volume 
increase  in  the  brain)  would  likely  lead  to  unexpected  and  unwanted  pressure 
fluctuations.  Another  hypothesis  is  that  the  boundary  conditions  applied  to  the  skull 
distort  the  shape  of  the  skull  and  cause  it  to  expand  slightly;  this  phenomenon  is 
similar  to  what  we  are  currently  investigating  in  the  Lagrangian  simulation,  which 
is  subject  to  comparable  loading  conditions. 

3.2  Effects  of  Mesh  Resolution  on  Results 

Within  discretization  methods,  it  was  often  the  case  that  the  resolution  of  the  mesh 
greatly  affected  the  magnitudes  and  distributions  of  the  results  given  by  our  models. 
An  effect  of  mesh  resolution  on  the  Lagrangian  simulations  was  that  the  magnitudes 
of  the  pressures  observed  varied  significantly  with  resolution.  Peak  pressures  at 
some  points  increased  as  resolution  increased,  and  others  did  not.  Figure  11  plots 
the  pressures  observed  at  3  of  the  tracer  locations  at  3  different  mesh  resolutions 
over  time. 
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Fig.  11  Pressure  over  time  by  mesh  resolution  at  3  tracer  locations,  Lagrangian  model.  The 
graph  legends  indicate  the  number  of  elements  used  in  each  simulation.  The  red  dashed  lines 
display  the  inflection  points  in  the  simulation — that  is,  the  approximate  points  at  which  the 
velocity  changed  directions.  Data  points  were  collected  every  1.0  ms. 
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The  peak  pressure  patterns  and  magnitudes,  although  not  entirely  congruent,  are 
generally  consistent  across  mesh  resolutions  and  show  a  trend  toward  a  lower- 
pressure  magnitude  than  is  given  by  our  simulations.  In  all  cases,  as  one  nears  the 
center  of  the  brain,  the  order  of  magnitude  of  the  pressure  values  decreases  slightly, 
as  would  be  expected  for  a  rotating  body.  The  remarkable  uniformity  in  shape  of 
the  curves  shows  that  not  all  the  pressure  effects  are  mesh  related  and  that  our  model 
is  capturing  the  physics  of  rotation  of  our  head  model  without  much  interference. 
The  resultant  pressures  in  the  Lagrangian  rotational  models  are  all  an  order  of 
magnitude  higher  than  the  corresponding  pressures  in  their  Eulerian  and  embedded 
counterparts.  As  mentioned  previously,  the  steady  increase  in  pressures  is 
completely  dependent  on  changes  in  the  volume  of  the  brain  that  may  in  turn  result 
from  errors  in  our  rotational  boundary  conditions;  investigation  into  the  cause  of 
the  volume  change  is  ongoing. 

In  the  Eulerian  simulations,  the  peak  pressure  values  were  similar  to  each  other,  but 
the  times  at  which  peak  values  were  reached  differed.  Also,  at  locations  closer  to 
the  skull,  the  pressure  patterns  observed  were  inconsistent  across  resolutions. 
However,  like  the  Lagrangian  simulations,  the  pressure  curves  appear  to  be 
converging  to  a  single  result  as  mesh  resolution  increases.  Figure  12  displays  this 
phenomenon. 
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Pressure  in  Brain  5  mm  in  from  Right  Inner  Edge  of  Skull,  Euler,  by 
Number  of  Mesh  Elements 
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Pressure  in  Brain  20  mm  in  from  Right  Inner  Edge  of  Skull,  Euler, 
by  Number  of  Mesh  Elements 
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Pressure  in  Brain  50  mm  in  from  Right  Inner  Edge  of  Skull,  Euler, 
by  Number  of  Mesh  Elements 
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Fig.  12  Pressure  over  time  by  mesh  resolution  at  3  tracer  locations,  Eulerian  model.  As  in 
Fig.  11,  the  graph  legends  indicate  the  number  of  elements  used  in  each  simulation,  and  the 
red  dashed  lines  display  the  inflection  points  in  the  simulation — that  is,  the  approximate  points 
at  which  the  velocity  changed  directions.  Data  points  were  collected  every  1.0  ms. 
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For  the  Eulerian  simulations,  the  primary  effect  of  mesh  resolution  appears  to  be 
the  phase  of  the  pressure  curve  as  opposed  to  pressure  magnitudes,  and  the  phases 
become  more  in  agreement  among  the  models  as  one  nears  the  center  of  the  brain. 
Since  the  Eulerian  method  averages  material  properties  over  mixed-material 
elements  (in  our  model,  along  the  skull/CSF/brain  boundary,  which  the  5-mm  tracer 
was  very  close  to),  mesh  resolution  would  be  expected  to  be  a  large  factor  in  the 
results  along  boundaries  between  materials  and  virtually  a  nonfactor  far  from 
mixed-material  elements.  As  mesh  resolution  increases,  the  area  over  which 
averaging  occurs  decreases,  potentially  improving  the  accuracy  of  the  Eulerian 
model.  One  of  the  most  important  pieces  of  Fig.  12  is  that  the  boundary  effects 
propagate  toward  the  center  of  the  brain;  across  resolutions,  the  highest  pressures 
were  located  20  mm  from  the  skull  as  opposed  to  the  normally  expected  5  mm  from 
the  edge  of  the  skull. 

In  creating  Figs.  11  and  12,  sets  of  results  for  5  different  mesh  resolutions  were 
initially  plotted,  but  this  number  was  reduced  to  3  to  allow  for  better  visualization 
of  results.  The  high,  low,  and  middle  resolutions  were  kept  in  the  graphs. 

The  unfortunate  effects  of  under-resolution  are  manifested  in  the  comparison  of 
embedded  grid  results  across  mesh  resolutions  as  shown  in  Fig.  13.  Note  that  the 
zero  pressure  for  these  simulations  is  atmospheric  pressure  instead  of  0.0  MPa. 
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Pressure  in  Brain,  5  mnn  in  from  Right  Edge  of 
Skull,  Embed,  by  Number  of  Mesh  Elements 
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Fig.  13  Pressure  over  time  by  mesh  resolution  at  3  tracer  locations,  embedded  model.  As  in 
Figs.  11  and  12,  the  graph  legends  indicate  the  number  of  elements  used  in  each  simulation. 
Data  points  were  collected  every  1.0  ms. 
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The  steepness  of  the  spikes  and  drops  in  this  figure  seen  in  the  meshes  with  lower 
resolutions  can  be  attributed  mainly  to  the  checkerboarding  in  the  simulation  that 
resulted  from  the  low  mesh  resolution.  The  amplitudes  of  the  peaks  decrease  with 
increasing  mesh  resolution,  implying  that  the  simulations  are  converging  to  a 
gentler  pressure  curve  that  hovers  very  near  atmospheric  pressure.  Like  the 
Lagrangian  simulations,  the  highest  pressures  seen  in  the  embedded  grid 
simulations  are  at  the  very  outside  of  the  brain. 

3.3  Efficiency  Analysis 

The  efficiencies  of  each  discretization  type  with  respect  to  simulation  size  are 
difficult  to  compare  across  all  3  methods;  however,  some  analysis  can  be 
performed.  Figure  14  displays  the  wall  times  of  the  Lagrangian  and  Eulerian 
simulations  with  respect  to  mesh  resolution.  All  simulations  were  run  on  Excalibur 
on  64  processors. 


Fig.  14  Wall  time  by  mesh  resolution  of  Lagrangian  and  Eulerian  methods,  rotational 
models 

In  finite  element  analysis,  it  is  usually  expected  that  Lagrangian  simulations  run 
faster  than  Eulerian  simulations  because  Lagrangian  models  do  not  have  to  account 
for  advection  across  element  boundaries.  This  is  definitely  the  case  for  our  models, 
and  it  appears  that  both  methods  provide  a  nearly  linear  run  time  with  respect  to 
input  (mesh)  size.  Because  of  the  time  constraints  of  our  project,  we  were  only  able 
to  obtain  results  at  3  different  mesh  resolutions  in  the  Lagrangian  case. 
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Unfortunately,  because  of  the  constraints  of  our  simulations  and  the  necessary 
restarts  that  had  to  be  performed,  the  exact  recorded  wall  times  of  the  embedded 
simulations  were  lost.  However,  we  can  approximate  that  the  embedded  simulation 
with  111,072  elements  ran  for  about  46  h  to  reach  40.0  ms;  the  Lagrangian 
simulation  with  approximately  twice  that  many  elements  ran  for  around  14  h  to 
reach  520.0  ms.  The  inefficiency  posed  by  the  embedded  model  was  a  major 
limiting  factor  on  the  quality  of  results  available  for  this  study. 

4.  Blast  Loading  Model  Analysis 


4.1  Results  across  Discretization  Methods 

Figure  15  depicts  the  results  of  the  blast  loading  simulations  using  each 
discretization  method  at  the  highest  resolution  we  were  able  to  run. 
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Pressure  in  Brain,  5  mm  in  from  Right  Edge  of  Skull,  by 
Discretization  Method,  Highest  Resolution 
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Fig.  15  Pressure  over  time  by  discretization  method  at  6  tracer  locations,  highest  resolution. 
The  number  of  elements  in  the  Eulerian,  embedded,  and  Lagrangian  simulations  are 
2,531,250;  2,532,018;  and  311,994;  respectively.  Data  points  were  collected  every  1.0  ms. 
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Pressure  in  Brain,  50  mm  in  from  Right  Edge  of  Skull, 
by  Discretization  Method,  Highest  Resolution 
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Fig.  15  Pressure  over  time  by  discretization  method  at  6  tracer  locations,  highest  resolution. 
The  number  of  elements  in  the  Eulerian,  embedded,  and  Lagrangian  simulations  are 
2,531,250,  2,532,018,  and  311,994,  respectively.  Data  points  were  collected  every  1.0  ms 
(continued). 
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Pressure  in  Brain,  20  mm  in  from  Top  Edge  of  Skull,  by 
Discretization  Method,  Highest  Resolution 
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Fig.  15  Pressure  over  time  by  discretization  method  at  6  tracer  locations,  highest  resolution. 
The  number  of  elements  in  the  Eulerian,  embedded,  and  Lagrangian  simulations  are 
2,531,250,  2,532,018,  and  311,994,  respectively.  Data  points  were  collected  every  1.0  ms 
(continued). 
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Unlike  in  the  rotation  problem,  the  magnitudes  and  distributions  of  pressures  in  the 
blast  loading  simulations  across  discretization  schemes  were  remarkably  similar. 
Each  simulation  clearly  showed  5  distinct  pressure  wave  reflections  within  the  skull 
over  the  course  of  the  simulation;  these  reflections  were  observed  at  almost 
identically  the  same  times  and  in  the  same  places  in  each  case.  The  magnitudes  of 
pressures  never  deviated  from  atmospheric  pressure  by  more  than  0.1  MPa  in  any 
of  the  cases.  The  primary  differences  among  the  sets  of  results  lie  at  the  points 
closest  to  the  brain/skull  interface,  as  expected.  The  embedded  grid  methodology 
displays  the  shallowest  wave  peaks  along  the  edge  of  the  brain;  the  air  required 
behind  the  skull  for  the  embedded  model  may  have  had  the  effect  of  dampening  the 
pressure  wave.  The  sharpest  and  tallest  wave  peaks,  as  well  as  clearest  pressure 
modes,  are  seen  in  the  Lagrangian  simulations;  since  Lagrangian  meshing  does  not 
have  mixed  material  zones,  the  averaging  of  properties  seen  in  the  other  2 
methodologies  may  be  dampening  results  by  providing  a  gradient  of  densities  over 
which  the  pressure  wave  has  to  pass  before  reflecting. 

4.2  Effects  of  Mesh  Resolution  on  Results 

Mesh  resolution  of  the  blast  loading  meshes  was  not  as  large  a  factor  in  the  pressure 
values  recorded  over  time  as  it  was  for  the  rotational  simulations.  Figure  16  plots 
the  results  for  the  Eulerian  blast  loading  model. 
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Pressure  in  Brain,  5  mm  from  Right  Inner  Edge  of 
Skull,  Euler,  by  Number  of  Mesh  Elements 
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Fig.  16  Pressure  over  time  by  mesh  resolution  at  3  tracer  locations,  Eulerian  blast  model. 
The  graph  legends  indicate  the  number  of  elements  used  in  each  simulation.  Data  points  were 
collected  every  1.0  ms. 
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Pressure  magnitude  and  wave  reflection  timing  are  extremely  consistent  across 
mesh  resolutions  for  the  Eulerian  simulations.  The  “shakiness”  observed  in  the 
values  across  resolutions  at  the  5-mm  tracer  particle  is  likely  due  to  the  mixed- 
material  elements  at  the  brain/CSF/skull  interface;  finer  resolution  of  the  CSF 
results  in  mixed  elements  comprising  a  smaller  portion  of  the  mesh,  which  should 
improve  the  accuracy  of  the  simulation. 

Mesh  resolution  appears  to  be  slightly  more  of  a  factor  in  the  embedded 
simulations’  results,  as  shown  in  Fig.  17. 
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Pressure  in  Brain,  5  mm  from  Right  Inner  Edge  of 
Skull,  Embedded,  by  Number  of  Mesh  Elements 
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Fig.  17  Pressure  over  time  by  mesh  resolution  at  3  tracer  locations,  embedded  blast  model. 
The  graph  legends  indicate  the  number  of  elements  used  in  each  simulation.  Data  points  were 
collected  every  1.0  ms. 
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However,  the  convergence  among  the  models  is  still  very  tight.  Like  the  Eulerian 
blast  simulation,  the  greatest  divergence  in  pressure  values  is  seen  closest  to  the 
skull.  In  the  case  of  the  embedded  simulation,  this  may  be  due  to  accumulated  error 
in  the  FEusion  calculation  that  is  performed  at  the  brain/skull  boundary. 

The  Lagrangian  blast  loading  models  have  the  greatest  results  discrepancies  with 
differing  mesh  resolution  (Fig.  18);  part  of  this  may  be  due  to  the  smaller  numbers 
of  overall  elements  than  their  Eulerian  and  embedded  grid  counterparts.  The  results 
appear  to  converge  over  time  toward  a  clearly  defined  pattern  of  5  wave  reflections 
in  the  head  over  the  course  of  the  run.  Once  again,  the  highest  discrepancies  are 
closest  to  the  skull.  There  are  no  mixed-material  zones  in  the  Lagrangian  models; 
however,  embedded  grid  methods  were  used  to  implement  the  background  air, 
some  of  which  lies  behind  the  skull.  It  may  be  that  the  discretization  of  this 
calculation  is  a  large  factor  in  results.  Another  possible  method  of  implementing 
the  background  air  is  to  use  a  conformal  mesh  surrounding  the  head  with  a  slide 
surface  between  the  air  and  the  skull.  Unfortunately,  there  was  not  enough  time 
available  to  create  a  conformal  mesh  during  this  study.  Further  investigation  into 
different  methods  of  implementing  the  background  air  would  be  required  to 
understand  this  phenomenon  more  fully. 
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Pressure  in  Brain,  5  mm  in  from  Right  Edge  of 
Skull,  Lagrange,  by  Number  of  Mesh  Elements 
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Fig.  18  Pressure  over  time  by  mesh  resolution  at  3  tracer  locations,  Lagrangian  blast  model. 
The  graph  legends  indicate  the  number  of  elements  used  in  each  simulation.  Data  points  were 
collected  every  1.0  ms. 
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4.3  Efficiency  Analysis 


An  interesting  pieee  of  the  blast  loading  simulations  is  that  contrary  to  what  is 
normally  expected  in  modeling,  the  Lagrangian  simulation  appears  to  be  the  least 
efficient  for  the  number  of  mesh  elements  used  (Fig.  19).  This  is  probably  due  to 
the  use  of  embedded  grid  methods  to  provide  the  background  air  for  the  run;  it  is 
recommended  that  a  methodology  that  does  not  use  an  embedded  methodology  be 
implemented  for  the  Lagrangian  run.  However,  the  Lagrangian  simulations  only 
required  one  embedded  interface  (between  the  skull  and  air),  while  the  embedded 
simulations  required  2  embedded  interfaces  (between  the  skull  and  air  and  between 
the  brain  and  skull).  In  this  planar  situation,  the  Eulerian  model  is  by  far  the  quickest 
run. 


Wall  Time  per  Cycle  by  Number  of  Mesh 
Elements,  Blast  Loading  Simulations 
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Fig.  19  Wall  time  by  mesh  resolution  of  each  discretization  method,  blast  models.  For 
clarity,  a  projected  linear  trend  curve  is  included  for  the  Lagrangian  model. 

5.  Conclusions  and  Recommendations 


We  have  seen  that  in  under  both  rotational  and  blast  loadings,  brain  trauma 
simulation  results  are  significantly  affected  by  the  discretization  method  used  in 
creating  the  mesh.  In  the  rotational  problem,  results  using  the  same  boundary 
conditions  were  disparate;  there  were  little  similarities  among  the  pressure  wave 
patterns  observed  in  each  case.  The  embedded  grid  model  appears  to  have  the  most 
spatially  similar  pressure  patterns  to  the  data  from  human  trials.'  However, 
although  the  preliminary  results  of  the  method  look  promising,  the  pressure  loss 
during  the  simulation  and  the  extremely  long  run  times  of  the  model  make  it  less 
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than  desirable  for  the  application  we  have  presented.  Investigation  into  the  cause  of 
the  pressure  loss  and  into  techniques  that  improve  the  efficiency  of  the  simulation 
are  of  paramount  importance  for  future  use  of  the  embedded  method  in  cases  when 
the  foreground  and  background  meshes  slide  against  each  other.  Also,  the 
embedded  model’s  usefulness  was  limited  by  the  low  mesh  resolutions  required  to 
complete  the  simulations  on  time;  for  further  investigation,  is  it  essential  that  a  finer 
embedded  mesh  be  used.  The  Lagrangian  rotational  simulation  requires  further 
study  as  well  to  determine  the  cause  of  the  volume  change  that  results  in  a  pressure 
buildup  that  deviates  from  expected  results.  Our  hypothesis  is  that  the  volume 
change  results  from  error  in  the  boundary  conditions,  but  this  hypothesis  requires 
verification. 

The  most  pressing  topic  for  further  investigation  across  all  methods  is  the  effect  of 
brain  structure  on  the  propagations  of  pressures,  stresses,  and  strains.  It  is  expected 
that  a  geometry  which  accurately  represents  the  folds  and  lobes  of  the  brain  will 
have  enough  internal  sliding  mechanisms  that  the  shear  stresses  will  remain  highest 
on  the  outside  of  the  brain,  which  is  very  different  from  our  results.  For  the  blast 
loading  problem,  all  3  methods  perform  similarly  and  efficiently  in  areas  far  from 
the  skull.  An  appropriate  next  step  would  be  to  conduct  a  study  that  provides  a 
clearer  physical  explanation  for  the  discrepancies  in  results  near  the  brain/skull 
interface.  Also,  methods  to  implement  the  Lagrangian  model  that  do  not  use 
embedded  techniques  to  implement  the  air  are  of  interest  in  future  work  both  to 
isolate  the  Lagrangian  method  within  the  simulation  and  to  improve  the  model’s 
run  time.  Unlike  the  rotational  problem,  it  is  not  anticipated  that  incorporating  an 
accurate  brain  geometry  will  affect  results  since  no  shear  forces  are  involved  in  the 
simulation. 
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Appendix.  Material  Parameters 


This  appendix  appears  in  its  original  form,  without  editorial  ehange. 
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Void 


voidinput  ssO  1000 


Air 

matinput 

czero  2.0  qfb  0.0  crq  0.1  pmin  -l.Oe-01 
eosvmin  l.Oe-03  eosvmax  l.Oe+03  vhlimit  l.Oe+05 
rho  1.18e-6  cvav  {  (720 . ) ^  (1 . 18e-6)  }  pOeosvmax  1 

evfrom  tp 

pO  1.01325e-l 
to  298. 

eO  2.533125e-01 
vO  1.0 
koinput 

isol  0  iform  5 
coef  0.4 


Brain 

matinput 

rho  1.04e-3  eO  0.0  vO  1.0 

pmin  -2.0e2  epsfail  25 

czero  2.0  crq  .1  qfb  0.15  linq  1 

cvav  2.926  eosvmin  0.3  eosvmax  2.0 

to  300.0 


msinput 
ysmodel  146 

shr_mod  13.e-3  wO  0.3  order  2  adv_log  1 
elasmodel  99 
hardmodel  299 
eosmodel  304 

rhoc2  2200  si  2.56  s2  -1.986  s3  .2268  gO  0.5  a  0. 
eosvmin  0 . 3  eosvmax  2 . 0 
ecO  -877.8 
emO  800. 
failmodel  499 


Cerebrospinal  Fluid 

#  yO  from  Bloomfield  et.  al . ^  1998  ->  assumed  viscosity  and 

density  of  water 

matinput 

rho  1.04e-3  eO  0.0  vO  1.0 

pmin  -2 . 0e2  epsfail  25 

czero  2.0  crq  .1  qfb  0.15  linq  1 

cvav  2.926  eosvmin  0.3  eosvmax  2.0 

to  300.0 
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msinput 
ysmodel  100 
elasmodel  0 
shr_mod  2660 
hardmodel  201 

yO  8.9e-7  ybet  0  n  0.1  ymax  5  rate_a  0  rate_b  1 
rate_exp  1 
eosmodel  304 

rhoc2  2200  si  2.56  s2  -1.986  s3  .2268  gO  0.5  a  0. 
eosvmin  0 . 3  eosvmax  2 . 0 
ecO  -877.8 
emO  800. 
failmodel  499 

Skull 

matinput 

rho  1.412e-3  eO  0.0  vO  1.0 
cvav  3.772 

msinput 
ysmodel  100 
elasmodel  0 

shr_mod  2660 
hardmodel  201 

yO  4  ybet  0  n  0 . 1  ymax  5 
eosmodel  304 

rhoc2  3891  si  0.94  s2  0  s3  0  gO  1.0  a  0. 
ecO  -1131.6 
emO  800. 
failmodel  499 
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